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We consider the problem of denning quantum integrability in systems with finite number of 
energy levels starting from commuting matrices. We argue that if the matrices depend on a (real) 
parameter, one can provide a definition from this feature alone, leading to specific results such as 
exact solvability, Poissonian energy level statistics and to level crossings. 



I. INTRODUCTION 

The field of quantum integrable systems, originally a 
somewhat abstruse topic, has received a great deal of at- 
tention in recent years with the realization that many 
typical models are realizable in atomic, mesoscopic and 
macroscopic many-body systems. In describing molecu- 
lar systems such as benzene, we deal with Hamiltonians 
defined on a finite-dimensional Hilbert space. These arise 
from a lattice of finite size, with spins, fermions or bosons 
populating the sites. In condensed matter systems too, 
one often studies finite size systems as a prelude to taking 
the thermodynamic limit or to obtain finite size correc- 
tions. Our concern in this work is with the characteri- 
zation and definition of quantum integrability in a finite 
dimensional context, a task that is considerably more 
delicate than the corresponding classical case. 

We may view typical condensed matter problems as 
models with a finite number of discrete single-particle en- 
ergy levels, and sometimes (but not always) with a fixed 
number of particles. The Hamiltonian is then 'just' an 
N x N Hermitian matrix, with a suitable N. This de- 
scription is, by design, far removed from its parentage in 
the space of many body models. Given such a Hamilto- 
nian matrix, can we say whether it is integrable or not? 
What is the precise notion of quantum integrability in 
this case; can we separate Nx N Hermitian matrices into 
two clearly distinct classes - integrable and otherwise? If 
yes, what are the consequences of quantum integrability 
that can be derived from its definition, i.e. the charac- 
teristic properties of such 'integrable' matrices? 

In this paper, we present a coherent view point that 
emerges from our recent exploration of such questions. 
We present a summary of our main results and also list 
some open questions that remain unanswered. To put 
this enquiry into context, note that in the thermody- 
namic limit N — > oo, impressive alternative approaches 
are available. For example Baxter's work on the spectra 
of infinite dimensional corner transfer matrices^ culmi- 
nates in the realization of extra symmetries and struc- 
ture arising in that limit, formalized by the Yangian 
approaclP. Another viewpoint is of geometric origin due 
to Sutherland 3 and emphasizes non-diffractive scattering 
as a true hallmark of quantum integrability. Much as we 
admire these powerful viewpoints, our chosen task of un- 



derstanding finite dimensional matrix systems takes us 
in quite a different direction. 

In classical mechanics a system with n degrees of free- 
dom is said to be integrable if it has n functionally inde- 
pendent integrals of motion that Poisson commute^, i.e. 
are in involution. An exact solution of equations of mo- 
tion follows from this definition, and it further can be 
shown that the motion is confined to invariant tori cut 
out in the phase space by the conservation laws. Unfor- 
tunately, a straightforward import of this elegant notion 
of integrability into quantum mechanics is problematic. 

The difficulties are at least two fold: firstly what is a 
degree of freedom? In quantum models with a fixed num- 
ber of particles, this number seems closest to the classical 
notion of a degree of freedom. However, in addition to 
the particle number there are other integer parameters 
to take into account such as the number of sites or the 
magnitude of spins. For example, in the Hubbard model 
with L sites do we count the number of electrons n e or 
the number of holes 2L — n e l Or does the Heisenberg 
model for L spin-| particles have the same number of 
degrees of freedom as that for L spin-5 particles? 

The second difficulty is in defining what constitutes a 
nontrivial integral of motion. From an elementary theo- 
rem in algebrgP, we know that an arbitrary N x N Her- 
mitian matrix commutes with N other Hermitian matri- 
ces that may be chosen as the projection operators along 
each component in the diagonal representation. We could 
equivalently express the matrix and all its commuting 
partners as a power series in some chosen non degener- 
ate matrix. A similar statement can be proved for general 
commuting Hermitian operators with discrete or contin- 
uous spectra™^. If now we are given a special Hamilto- 
nian matrix, and told that it is quantum integrable (in 
some specific sense that is yet to be defined), possess- 
ing n equally special commuting partner matrices; we 
might (rather stubbornly) ignore this extra information 
and compute the N elementary commuting partners of 
the given matrix as per the earlier prescription. Don't we 
then have too many partners? What, if any, is the dis- 
tinction between these two sets of commuting partners? 
From such questions and considerations one seems to be 
forced to call integrable either all or none Hamiltonians 
with finite Hilbert spaces. 

Our goal here is to attempt a resolution of these dif- 
ficulties and provide a practically useful, simple and at 
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the same time rigorous definition of quantum integrabil- 
ity for Hamiltonian matrices. Our interest is of course 
not just in the definition per se but in a clear delineation 
and derivation of properties of quantum integrable sys- 
tems from the definition. Characteristics normally at- 
tributed to integrable models include: exact solution for 
their spectra, Poisson level statistics, crossings of levels 
of the same symmetry in parameter-dependent Hamil- 
tonians etc. However, they are never derived from one 
another or within some unified framework, but each prop- 
erty has to be established independently and on a model 
by model basis. 

Consider, for example, the Id Hubbard model. Shortly 
after the model was exactly solved via Bethe's Ansatz^, 
a separate rigorous studylEI of its parameter independent 
symmetries revealed numerous level crossings violating 
the Wigner von Neuman non-crossing rule for eigenval- 
ues in generic systems. These findings can be rationalized 
by the notion that copious level crossings are a hallmark 
of quantum integrabilty, and hint at parameter depen- 
dent conservation laws in such a model. Much later, one 
of u d 11 * 12 ! found an embedding of the Id Hubbard model 
into the Yang Baxter scheme, thereby displaying an ex- 
plicit set of parameter dependent conservation laws. Not 
surprisingly, a numerical analysis^ of level spacings finds 
Poisson statistics that crosses over to GOE when 'inte- 
grability' is destroyed by adding an off-site Coulomb re- 
pulsion. Further worlfHlby us on making a causal link be- 
tween the conservation laws and the level crossings in the 
Id Hubbard model has been a fruitful source of insights 
that we have subsequently explored and report here. 

Often, one of the above characteristics is singled out 
and adopted as a definition^. It seems to us that while 
level statistics and crossings are useful, and even powerful 
tests of integrability, these are hardly suitable as defini- 
tions. These criteria arguably encompass a broader class 
of systems than those normally thought integrable^. It 
is also difficult to see how, starting from either of them, 
one could derive other properties, e.g. an exact solu- 
tion. Universal statistics emerges only in the limit of 
large matrices, N — > oo, or for an ensemble of matrices, 
while here we are looking for a notion that also works for 
a stand-alone Hamiltonian matrix with fixed N. In ad- 
dition, there are exceptional points in parameter space 
of systems usually recognized as integrable, where the 
statistics is non-Poissonian, e.g. for Gaudin mag netsP. 
Starting from a well-defined notion of integrability that 
we propose below, we indeed find that there are such de- 
viations in a much broader class of systems and explain 
their origin. Similarly, we find examples that fail the 
level crossing test - such as 'accidental' degeneracies in 
non-integrable systems and (rare) instances of integrable 
systems without crossings. 

It turns out that such a simple and yet well-defined 
notion can in fact b e formulated as we recognized in a 
series of paperJSE^Ell. rp^g ma j n i^ea, j g ^ consider the 

dependence of commuting operators on a real parame- 
ter, which we denote u. This is typically an interaction 



constant or an applied field, e.g. Coulomb interaction 
constant in the Hubbard model, magnetic field in Gaudin 
magnets^, pairing strength in the BCS and anisotropy in 
the XXZ Heisenberg models, etc. Besides usual space- 
time and internal space symmetries, which are parameter 
independent, these models have it-dependent conserva- 
tion laws (alternatively termed dynamical symmetries or 
conserved currents) as discussed above. For the BCS and 
Gaudin models the Hamiltonia n and all conserved cur- 
rents are linear in the paramete r^ 2 * 23 !. In the XXZ and 
Hubbard models the Hamiltonian and the first dynami- 
cal symmetry are linear i n u; t he rest are polynomials in 
it of order two and higheiP^IIEIj] 



II. DEFINITION 

In most of the examples of condensed matter models 
believed to be quantum integrable, there is a set of com- 
muting operators linear in u including the Hamiltonian. 
This seems to be the most common situation in param- 
eter dependent integrable models and leads ua^HSl] to 
consider operators of the form H(u) = T + uV , where 
T and V are N x N Hermitian matrices. A key ob- 
servation is that fixed u-dependence implies a natural 
well-defined notion of a nontrivial integral of motion. A 
typical (e.g. randomly generated) H(u) commutes only 
with (a + bu)l + cH(u), a trivial operator linear in u. 
Here a, b, and c are real numbers and 1 is the identity 
operator. 

The requirement that there exist a nontrivial commut- 
ing partner linear in u severely constrains the matrix ele- 
ments of H(u). In fact, as we will see below, a real sym- 
metric H(u) that has such a partner is fixed by less than 
(N - l)(JV + 8)/2 real parameters. In contrast, N(N+1) 
real parameters are necessary to specify the matrix el- 
ements of T and V for a generic H(u), indicating that 
matrices with fixed parameter-dependence split into two 
non-overlapping sets - those with nontrivial commuting 
partners and those without. 

Moreover, there is a natural classification of integrable 
families according to the number n of linearly indepen- 
dent commuting operators. Specifically, we define an in- 
tegrable family as a vector space of n linearly independent 
N x N Hermitian matrices 

H k {u) = T k + uV k such that, 
[H l (u),W{u)\ = for all it and i,j = l,...,n. (1) 

In addition, we impose an (optional) condition that 
H l (u) have no common u-independent symmetries - 
there is no constant matrix Q al) such that 

[fl,H l (u)] — for all u and i. If there are such sym- 
metries, H l (u) are simultaneously block-diagonal and 
Eq. ([I]) reduces to that for smaller matrices (blocks) with- 
out u-independent symmetries^. 

Linear independence means that ^2 i ciiH t (u) = (a + 
bu)l with real ai if and only if all = and a = b = 0. 
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Note that, for convenience, we chose to separate multiples 
of identity from our list of nontrivial commuting matrices 
H k (u). Finally, n is defined as the maximum number of 
nontrivial independent commuting matrices in a given 
family, i.e. any H (u) = T + uV that commutes with all 
H l (u) can be written as 



H(u) =J2 a kH k {u) + (a + bu)l. 



(2) 



Thus H l (u) act as basis vectors in the n-dimensional vec- 
tor space V n of commuting matrices, defined up to a mul- 
tiple of identity matrix 1. 

We propose the following classification of integrable 
matrices linear in a parameter into types. The max- 
imum possible number of linearly independent N x N 
commuting Hermitian matrices is n = N — 1, not count- 
ing multiples of the identity matrix. We call this a type 1 
or maximally commuting family of matrices. Similarly, 
families with n = N — 2 independent commuting H l (u) 
are termed type 2, and a general type M family is de- 
fined through n = N — M. The maximum value of M 
is M = N — 2 when there are two nontrivial members in 
the family in addition to (a + bu) 1. 

We argue that the following properties are direct con- 
sequences of the above definition of quantum integrabil- 
ity: (1) an exact solution for the eigenspectra of H l (u) in 
terms of roots of a single algebraic equation, (2) H l {u) 
satisfy Yang-Baxter equation, (3) eigenvalues of H l (u) 
typically (but not always) cross as functions of u\ the 
number of crossings depends on both N and M (A) eigen- 
values have Poisson statistics in N — > oo limit except for 
some special cases of certain measure zero in the space 
of all integrable families. Below we prove statements (1) 
through (3) for type 1 families (in this case crossings 
are always present) and comment on similar results of 21 
for other types. Poisson statistics will be demonstrated 
numerically in a separate paper; here we briefly discuss 
some of its main findings. 

Eq. ([I]) can be solved for matrix elements of H l {u) at 
least for some types of int egrab le families. All type 1 
families were constructed ir l 18 * 19 [ IrP^ all type 2, 3 and 
some type M for arbitrary M were obtained. To facilitate 
further discussion, let us first cast Eq. ([!]) into a different 
form. Using H k (u) = T k + uV k and equating to zero 
terms at all orders of u, we obtain 

[T\ V 3 ] = [T\V% [T\T J 'j = [V\ V J ] = 0. (3) 

It is convenient to choose the basis in the target Hilbert 
space to be the common eigenbasis of the mutually com- 
muting matrices V 1 . The first commutation relation in 



Eq. (|3j in this basis reads T| m /(4 - = T{J{£ k - 
d J m ) = A km , where d\ are the diagonal elements of V 1 . 
This implies that T % can be written as 

T = W*+ [V\A] (A) 

where W % is a diagonal matrix. Note that the an- 
tihermitian matrix A is the same for all members of 



the family and is therefore independent of the basis in 
V n . Now the commutation relations [V 1 , V J ] = and 
[T*, V 3 *] = [T J , V 1 ] are satisfied. The remaining equation, 
[T\TJ] = 0, takes the form 

[[V\A],[V\A]] = [\y^A],W*] - [[V\A],wq . (5) 

Before we proceed with type 1, let us discuss the num- 
ber of parameters involved in constructing a generic real 
symmetric integrable family. 

We can generate integrable families of various types 
by solving Eq. ^ numerically^. The algorithm is as 
follows. First, we arbitrarily specify some diagonal ma- 
trices V 1 , V j , W % and W j (4N real inputs) and solve 
for the antisymmetric matrix A. This yields a dis- 
crete set of solutions for A and two commuting matrices 
H id (u) = W 1 ^ + [V i,j ,A] + uV^ for each A. To deter- 
mine, the remaining basic matrices H k (u) in the family 
with a given A, we take j — k in Eq. (pM) (with V" 1 , W l 
and A obtained before) and solve this equation for V k 
and W k . 

Some choices of the 4N inputs in the above procedure 
produce the same commuting family. In fact, as seen 
from Eq. ([2 1, there is 2(n+2)-parametric freedom (values 
of etfc , a and b) in picking two matrices within the family. 
This means that we can fix 2(n + 2) out of AN parameters 
by taking linear combinations within the family, i.e. by 
going to a different basis in the vector space of commuting 
matrices. In addition, Eqs. ^ and Q are invariant with 
respect to rescaling V 1 — > aV l and A— > A/ a, which fixes 
one more parameter in V % . This leaves 4./V— 2(?i+2) — 1 = 
2N + 2M — 5 real parameters to specify a generic type 
M family. 

To further select a particular matrix H(u) = T + uV 
within the family, one needs to pick n + 2 coefficients in 
Eq. |2| , in addition to 2 N + M — 5 that specify the com- 
muting family, i.e. a general type M matrix involves 
AN — n — 3 = 3./V + M — 3 real parameters. Since 
M < N — 2 the maximum number is AN — 5. This 
counting has been done in the common eigenbasis of V 1 . 
Going to an arbitrary basis in the target Hilbert space 
adds another N(N — l)/2 real parameters for a general 
orthogonal transformation. We see that the total number 
of parameters is less than (N + 8)(N — l)/2, about half 
of N(N + 1) for a non-integrable real-symmetric matrix 
of the form B + uC. 



III. TYPE 1 

The 'master' equation ([5| is simplest for type 1. In this 
case there are N linearly independent H l (u) (including 
the identity) and consequently iV linearly independent 
diagonal matrices V 1 . By taking linear combinations we 
can go to a basis in the vector space Vn-i such that 
= 4 = S ik . Eq. (g) reads f* k = W] 3 - W\ k = 
—AijAik/Ajk- These equations are consistent when /j fe + 
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fli + fij = 0' yielding the following four index relational 
AijAjkAkiAu + AikAkiAijAji + A a AijA ]k A ki = 0. (6) 
The most general solution of this equation i d 19 * 20 * 



A jk = 



(7) 



where real Si and complex 7^ are unrestricted parameters 
that fix the commuting family. 

Next, using Eq. (JsJ) and = 5 k jSik, we determine 
W l from Eq. ([s]), which is linear in W l , and T l from 
Eq. Q. The most general member of a type 1 family, 
H{u) — Yli=i diH l (u) with arbitrary real d^, is 



Wl mn = 7m7rt ( ~ — ) , m^n, 



[#( U )Lm = U 



rfm - 51 IT* 



2 ( dm &k 

£■171 ~~ 



(8) 



Note that d m are eigenvalues of V by design. In particu- 
lar, for basic operators H l (u) we have dk = Sik, i.e. their 
nonzero matrix elements are 



[**(«)]« = 



Mi 



Sij — Sik 



(9) 



A similar construction is possible for real symmetric in- 
tegrable families of arbitrary type A$^. It involves M + 1 
new parameters in addition to 7,, £j and dj, which play a 
role similar to that in Eq. ^ . However, this construction 
produces all real symmetric families only for M = 1, 2, 3, 
while for M > 3 it yields only a subset of such families. 

A correspondence between real symmetric type M 
families and compact Riemann surfaces of genus g > 
M— 1 was conjectured irPS The construction oW^ {ansatz 
type M) produces families that correspond to hyperel- 
liptic Riemann surfaces of genus g — M — 1. All Rie- 
mann surfaces of genus 0, 1, and 2 are hyperelliptic 
which explains the completeness of the construction for 
M = 1,2,3. 



IV. EXACT SOLUTION 

The exact spectra of type 1 matrices were obtained 
irP^l. The components of an eigenvector (column) v m (u) 
of H(u) given by ^ are 



[vm{u)} 



7i 



3 A — £■ 



with respective eigenvalue 



JY 



E m (u) = 



k=l 



d k \lk\ 2 



(10) 



(11) 



where the A^, i = 1, 
algebraic equation 



f(K, 



. , N are determined from a single 



E 

i=i 



I Til 



(12) 



These equations can be verified directly by evaluating 
J2 3 [H(u)]ij[v m (u)]j. Ansatz type M (see the last para- 



graph in Sect. Ill ) families have an exact solution in terms 



of a single equation similar to Eq. ( 12 ) 



All N roots A, of Eq. ( 12 1 are real. Indeed, /(A) -> +00 



as A — >• ej, and /(A) — > —00 as A ^ £ k+i> wnere £ fe are 
ordered, E\ < £ 2 < • • ■ < £/v- It follows that /(A) = u 
has a real solution between Sk and £&+i for any k, i.e. 
£m < A m < £m+i- The last root Xn > £jv for it > and 
A at < £at + i = £1 for u < 0, where we identify indices 
that differ by a multiple of AT. 



V. LEVEL CROSSINGS 

First, we show that any type 1 matrix H(u) = T + 
has at least one level crossing^. To this end, let us 
analyze the evolution of eigenvalues E m (u) with u. We 
observe that X m — ¥ e m as u — > +00. In this limit the 



main contribution to Eqs. (11) and (12 1 comes from the 
j = m term, yielding E m — > \u\d m . Similarly, we obtain 
x m — > £ m +i and E m — > —\u\d m+ i for u — > —00. It is 
not surprising that eigenvalues of H(u) tend to ±|u|dj, 
since dk are eigenvalues of V and uV dominates H(u) for 
large u. What is important however is that we know to 
which particular udk a given E m (u) tends in both limits. 
Symbolically, we can write k — > k — 1 (mod N) meaning 
the eigenvalue goes from — \u\dk on the left (u — > —00) 
to |u|dfe_i on the right (u — > +00). 

The presence of levels crossings can now be proved by 
contradiction. Suppose there are no crossings. Since 
eigenvalues are continuous functions of u, this implies 
that their ordering must be the same at all u. The top 
level must connect the largest eigenvalue at u — » —00 
to the largest eigenvalue at u — > +00, the bottom level 
goes from the lowest eigenvalue at u — > — 00 to that at 
u — > +00, etc. Let dk be ordered as^di < dj < ■ ■ ■ < d m . 
Then, the largest (lowest) eigenvalue at u — > —00 is 
— \u\di (— \u\d m ) and the largest (lowest) eigenvalue at 
u — > +00 is \u\d m {\u\di), i.e. we have i — > m and m — > i 
for the top and bottom levels, respectively. On the other 
hand, according to the k — > k — 1 rule established above 
this implies m = i — 1 (mod N) and at the same time 
i = m — 1 (mod N). We obtain = 2 (mod N), which 
does not hold for any N > 3, i.e. the above assumption 
that levels do not cross cannot be true. Thus, at least 
one level crossing is inevitable. 

Allowed values of the total number of crossings n x in a 
type 1 matrix can be determined by analyzing the repre- 
sentation of an arbitrary type 1 matrix H(u) in terms of 
powers of any other nontrivial H(u) that belongs to the 
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FIG. 1: A schematic energy level diagram for N — 5 demon- 
strating maximum number of crossings (TV — 1)(N — 2)/2. The 
diagram is for the ordering di > di > ■ ■ ■ > djv and uses the 
rule k — > k — 1 derived in the text. 



same commuting family. One finds^ that this expansion 
is necessarily of the form 



JV-l 



H(u) = £ 



Q m (u) 
P H (u) 



[H(u)]' 



(13) 



where P H (u) is a polynomial in u of degree (N — 1)(JV — 
2)/2 with real coefficients that depend on matrix ele- 
ments of H(u) only. Q m (u) are polynomials in u of order 
m — 1 lower than P (u). This expansion breaks down at 
a given value of u — Ux only when H{u) has a crossing 
at u = u x and P H (u x ) — 0. Thus, crossings of H(u) 
occur at the roots of P H (u) (see^Sfor a detailed proof). 
We conclude that the maximum number of crossings in 
a type 1 matrix is 



(N-l)(N-2) 



(14) 



This upper bound is realized e.g. for matrices ([8]) such 
that d\ > g?2 > • • • > d,N, see Fig. [T] n x has a definite 
parity, that of n™ ax , because the coefficients of P H (u) are 
real and its complex roots therefore come in conjugate 
pairs. For example, 4x4 type 1 matrices have either 1 
or 3 crossings, n x — 2, 4, 6 for 5 x 5 etc. 

Crossings are also ubiquitous in type M > 1 in- 
tegrable families. It is possible however to deliber- 
ately engineer 4x4 type 2 matrices without crossings^. 
Other than that, we do not have rigorous results for 
higher types. Empirically, one finds 2 ^ that the expan- 
sion ( |13[ ) still holds, but the degree of P H (u) and there- 
fore the maximum number of crossings is reduced to 
(N - 1)(N - 2)/2 - g, where g > M - 1 is the genus of 
the corresponding Riemann surface, typically g = M — 1 
(see the end of Sect. III). 

Interestingly, one can also construct 4x4 real sym- 
metric matrices A + uB that have crossings but no u- 
independent symmetry and no commuting partners lin- 



ear in u (this is impossible for N = 3), se e! 14 * 19 ! f or de- 
tails. This indicates that either 'accidental' degeneracies 
are possible or one can meaningfully introduce nonlinear 
in u conserved currents in such cases. 



VI. YANG-BAXTER EQUATION FOR THE 
MATRIX MODEL 



We next show how the type 1 matrices can be fit into 
the Yang Baxter formulation of integrable systems, by 
displaying matrix objects S, Eq. (17 1, that play the role 



of scattering amplitudes. Thus we consider a purely ma- 
trix model in N dimensions with states \i), 1 < i < N, 
the projection operators 7r,-j and the identity matrix 1: 

1 = 22 Kg 

i 

With i j, let us define a 'dressed' permutation operator 

J^Li + 7l 7.: faj + Tift) - 



n 



7i 



l7j|V 



If we set 7i — > , II reduces to the permutation operator 

Hy -> fti, (15) 

where acts as: 

A#) = fob') + *jfc|t) + (i - fo)(i - Ml fc >- (is) 

Let us introduce j as a coupling parameter and also 
the composite parameter x = (e, 7) so that the scattering 
operator Sy is defined as: 



\ T .\ = <?.. - (£j - + 2g % 



(17) 



The action of a particular scattering matrix SVim on 
the relevant states is given by 



S nm \k) = \k), [k ^ n,m) 
S nm \n) = t+ m \n) + r nm \m) 
S nm \m) = t~ \m)+r nm In) 



(18) 



with the reflection (r) and transmission (t) amplitudes 
defined by 

257«7m 



nm ~ 


- £n " 


\-g{\jm\ 2 - 


f |7n 


2 Y 


-\- 


&n 


(-5(l7«| 2 - 


l7m 


2 ) 


nm 


~~ &n 


f"5(l7m| 2 - 


f l7n 


2 )' 


— £m 


- £n - 


Ml7m| 2 - 


" |7n 


2 ) 


nm 


- £n " 


^g{\im\ 2 - 




2 )' 



(19) 
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We state the two operator relations that are needed, 
and readily verified, for any three indices 1, 2, 3: 

ilia n 32 n 31 = n 31 n 32 n 12 , (20) 
[n 12 , (H31 + H32)] = 0. (21) 

With Cij as constants we can now verify the condition 

(c 12 i + g n 12 ) (c 32 i + 9 n 32 ) (C31I + 9 n 31 ) = 
(Caii + 9 n 31 ) (c 32 i + 9 n 32 ) (c 12 i + 9 n 12 ), (22) 

provided C s satisfy the triangle law: 

C 32 = C 3 i + Ci2- (23) 



Dividing by a suitable constant and consulting Eq. ( 18 ) 



we therefore verify the Yang Baxter equation for the S 

SikSjkSij = SijSjkSik, (24) 
and the initializing condition 



lim S l j(xj\x i ) 



(25) 



We consider the N sites (indices) and add an auxiliary 
index a (that is outside the original space of N states) 
so that Wi a , TT a i, ir aa are added to the list of operators. 
Now define a (monodromy) matrix 

T a = B a (x a ) Snoc Sn-i a - ■ ■ Sx a . (26) 

This is a function of all the indicated variables: 

T a = T a (x a \{xi,X2,.-.,x N }), (27) 

and also a boundary field term B a (x a ). For the boundary 
term (twist) to give commuting operators, this term must 
be chosen to satisfy the conditional 



[S a p,B a B ] = 0, 



whereby we choose 



9 



U 



(29) 



with the freedom of an arbitrary parameter u. 

The transfer matrix is obtained by tracing over the 
auxiliary index a 



T(x a \ {xi,x 2 , ■ • • ,x N }) = 22(a\T a \a) 



(30) 



Using Baxter's classic proofs for commutation of trans- 
fer matrices, slightl y gen eralized to the case of twisted 
boundary condition d 22 * 32 1 , we conclude that 



[T(x a ),T(x )]=O, 



(31) 



for arbitrary x a and xp while holding {xj} g, u fixed. 



We note that Eq. (31 1 is also valid as x a — > Xi and hence 



conclude [Tj,Tj] = 0, where 

Tj = ^lim T(x a \{xx,x 2 ,...,x N }) 

Bj (xj )S NJ {xj , x N ) . . . Sj+ij (xj ,x j+ i). (32) 



An expansion in powers of the interaction strength g pro- 
duces the currents: 
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T, = l + ^^( U ) + 0(^), 



(33) 



with 



Considering terms of order 0(g) in [Ti,Tj], we conclude 

[H i (u),H j (u)]=0. (34) 

Note that H l (u) are the exactly the basis operators of 
Eq. m written in terms of projection operators 2 ^!. 

VII. LINKS TO VARIOUS MODELS 

First, note that one can choose an arbitrary Hermitian 
matrix V and still have H(u) = T+uV to be a member of 
any given type 1 family. Indeed, Eq. ^ is written in the 
eigenbasis of V. Both this basis and the eigenvalues of V 
- parameters dk - can be chosen arbitrarily. By symme- 
try one can instead choose T at will, though this is not 
apparent from Eq. But as soon as e.g. V is fixed, T 
is severely constrained - one is left with only 2N param- 
eters 7i,£j to specify its matrix elements. This ability to 
choose either V or T arbitrarily means in particular that 
any it-independent Hamiltonian can be 'embedded' into 
a type 1 family in many different ways. For example, 
one can choose V to be the isotropic Heisenberg (XXX) 
model, or the Haldane-Shastry model, and find T so that 
T + uV belongs to a given type 1 family. 

Type 1 integrable families are closely related^ to 



(28) Gaudin magnets^ hi = Bs*+J2j^ti §iSj(ej— e^) 1 , where 
Sj are quantum spins of arbitrary length Sj . In the sector 
with (conserved) S z = • sf equal to its maximum (min- 
imum) possible value less (plus) one, hi are N commuting 
N x N matrices which form a a type 1 family with u = B 
and 7 2 = Sj. The BCS model is obtained^ as J^j £ jhj 
for 7 2 = Sj = 1/2 and a replacement u = B — s- 1/g, 
where g is the dimensionful BCS coupling constant. 

Some blocks of the ID Hubbard model characterized 
by a complete set of it-independent symmetry quantum 
numbers are type 1 matrices, though most blocks are type 
M > l 21 . A similar typology can be developed for e.g. 
the ID XXZ Hamiltonian and other sectors of Gaudin 
and BCS models using the method o^ for determining 
the type of parameter-dependent matrices. Interestingly, 
this implies that at least in some blocks there is an exact 
solution in terms of a single algebraic equation - Eq. ( 12 ) 
or a similar equation for higher types^H- a vast simplifi- 
cation as compared to Bethe's Ansatz. 

One can also construct fermionic (bosonic) 
Hamiltonians^Sl out of type 1 matrices as 



H = y^\H(u)] mn al n a n , 



(35) 
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where a n are the usual fermionic (bosonic) destruction 
operators. [H\,Pi2\ = as long as the corresponding 
matrices Hi(u) and H 2 (u) commute, i.e. belong to the 
same family. 



VIII. LEVEL STATISTICS 

We have performed an extensive numerical study of 
level statistics of type 1 and higher type matrices^ for 
various choices of parameters. Almost in all cases the 
statistics is Poissonian for JV 3> 1 with high accuracy, 
even when we deliberately attempt to adjust the param- 
eters to get a different statistics. Let us briefly describe 
the main results e.g. for the level-spacing distribution. 

It is convenient to redefine the parameter u —¥ 1/x 
and replace T + uV —> V + xT. To get a proper large 
JV limit one has to make sure that the scaling of param- 
eters dk, £fe ; 7fc in Eq. ([8| and x with JV is such that the 
eigenvalues of V and xT scale in the same way for large 
JV. As discussed above, the matrix V is arbitrary, so 
at x = one can have any admissible level statistics. 
Consider, for example, three representative cases: (a) 
V is a random real symmetric matrix with independent 
identically distributed matrix elements Vjk for j < k, 
(Jo) eigenvalues dk of V are independent uniformly dis- 
tributed random numbers, and (c) dk display level attrac- 
tion, P(s) = as" exp(-6s 1+ ") with -1 < u < 0. The 
level-spacing distribution P(s) for V is Wigner-Dyson 
P(s) = 2ase~ as in (a) and Poissonian P(s) = e~ s in 
(b). As soon as V is chosen, T is no longer arbitrary and 
we find that it has Poissonian P(s) in all three cases for 
all choices of parameters jk and £k we considered as long 
as £fc and dk are uncorrelated. 

Specifically, motivated in part by the BCS and Gaudin 
examples discussed in the previous section we took: (1) 
£fe that are also eigenvalues of a random matrix and 
7fe = const independent of k, (2) same as (1) but with 
random uncorrelated 7& and (3) independent uniformly 
distributed Ek and 7/- = const. For all these choices the 
level-spacing distribution for T is very well approximated 
by Poissonian P(s) = e~ s , where s is the level-spacing 
in units of the mean level-spacing. In case (a) above the 
level statistics of H(x) = V + xT at x = is Wigner- 
Dyson, but we find that it crosses over to Poisson at 
I or I ~ 1/JV and remains Poisson for larger \x\. Case (c) 
is analogous to (a) - a crossover to Poisson behavior at 
|a?| ~ 1/JV. In case (b) the statistics is Poissonian for all 
x. Similar behavior is found in spectral rigidity. We con- 
clude that one can arrange for any statistics at a given 
value of the parameter x = xq, but this becomes an iso- 
lated point in N — > 00 limit, while for x ^ xq integrability 
as defined in Sect. ITT1 enforces Poisson statistics. 

The only exception to Poisson statistics other than at 
an isolated value of x we were able to identify is when 
parameters dk and Sk are correlated, so that dk = f(£k), 
where /(e) is a smooth function of e in JV — ► 00 limit 35 
and 7^ = const. This is the case in e.g. the BCS model 



where dk = £k (see above). In such cases the statistics 
is distinctly non-Poissonian and, moreover, in case (a) 
above, for example, P(s) crosses over at x = 0(N°) 
( ) 1 1 from the Wigner-Dyson P(s) — 2ase~ as to a more 
repulsive distribution P(s) oc s 4 for small s. The repul- 
sion is softened by randomizing 7^. More importantly, 
the statistics quickly becomes Poissonian when the corre- 
lation between dk and £k is destroyed, dk = f(sk)0-+T)k), 
where r\k are random. We find Poisson distribution al- 
ready for r\k = 0(1/N) at x = 0(1), see alscP^ for a 
similar study of Gaudin model, dk = f(sk) define ex- 
ceptional 'surfaces' of certain measure zero in parameter 
space. This seems analogous to the harmonic oscillator 
exceptions to the Poisson distribution in classical inte- 
grable systems^. There too one finds increased level 
repulsion for oscillators instead of Poisson P(s). 

Some of these numerical observations can be under- 
stood using perturbation theory. Energies to the first 
order in x are given by the second equation in pi , where 
we set |7j| 2 = 1/JV to achieve proper scaling for large JV 
as discussed above^. We have 



E m (x) 



X ^ — "v 

1 * . , \ c r 



(36) 



The first term comes from V , which we take to have 
Wigner-Dyson P(s), the second - from T, which is deter- 
mined by the integrability condition ([3| and whose level 
statistics we do not control. Let us estimate x at which 



the two terms in Eq. (36) become comparable. Without 
loss of generality we can take dk = O(N ) = 0(1) and 
we must also take £k = 0(1) so that T and xV scale 
in the same way for large JV. Suppose £k are ordered 
as £\ < £2 < ■ ■ ■ < £n- When dk and £k are uncorre- 
lated d m — dj is O(l) when j is close to m, i.e. when 



(£ m — £j) = 0(1 /N). The second term in Eq. (36 1 is 
then a:c m lnJV, where c m — 0(1) is a random number 
only weakly correlated with d m . 

If we now order d m , c m in general will not be ordered, 
i.e. if d m+ i > d m is the closest level to d m and there- 
fore (d m +i — d m ) = 0(1 /JV), the corresponding difference 
(c m +i — c rn ) = 0(1). The contributions to level-spacings 
from the two terms in Eq. ( 36 1 become comparable for 
x = x c s» l/(JVlnJV). It makes sense that the second 
term introduces a trend towards Poisson distribution be- 
cause it is a (nonlinear) superposition of £k and dk - 
eigenvalues of two uncorrelated random matrices. Thus, 
we expect a crossover from Wigner-Dyson to Poisson dis- 
tribution at x = x c . 

This argument breaks down when dk = f(£k), since in 
this case (d m -dj) = 0( 1/JV) when (e m -Ej) = 0(1/ N). 
The two terms in Eq. ( 36 ) become comparable only at 
x = 0(1) in agreement with the numerics for this case. 
Moreover, the second term no longer trends towards Pois- 
son statistics. Relaxing the correlation between dk and 
Ek with dk = f(£k)(l + T)k) and going through the same 
argument, one expects a crossover to Poisson statistics 
&tx = 0(1) for rjk = 0(1/N). 



IX. DISCUSSION 

A distinct feature of the notion of quantum integrabil- 
ity proposed in this paper is the parameter (u) depen- 
dence. This is in contrast to the classical notion that 
does not require any such dependence. We however find 
it necessary to be able to quantize the classical definition 
in a meaningful way. In this sense our definition is more 
demanding than its classical counterpart. 

On the other hand, there is another distinction from 
the classical notion that makes our definition seem less 
stringent - even a single nontrivial integral of motion 
linear in u is sufficient to declare the model integrable. 
This was originally motivated by the absence of a well- 
defined analog of the number of degrees of freedom in 
quantum mechanics. There is evidence however that the 
presence of a single such integral is actually much more 
consequential in the quantum case. For example, for any 
H(u) = T + uV with a single nontrivial commuting part- 
ner I(u) = K + one also finds numerous (about N 
linearly independent) integrals of motion quadratic in u 
that commute with both H{u) and I(up^. We note that 
a generic matrix of the form A+uB has no quadratic inte- 
grals other than its own square and a multiple of identity. 
One of the implications of this is that, for example, cur- 
rents of higher order in the parameter in the ID Hubbard 
and XXZ models might follow from the linear ones, i.e. 
be in some sense trivial given the linear integral. These 
are however open questions that require further research. 

One can also consider e.g. a situation when a Hamil- 
tonian of the form H(u) = T + uV or T + uV + u 2 W 
has no nontrivial integrals linear in u, but a number of 
quadratic ones. Such systems exist and can also be classi- 
fied and, at least some of them, explicitly parametrized. 
Nevertheless, this 'higher order in the parameter' inte- 
grability seems less relevant as in most physical examples 
of parameter-dependent integrable models one is able to 
identify a parameter such that the Hamiltonian and at 
least one of the currents are linear in it. 

It is worth commenting on the relationship between our 
viewpoint and the usual set of 'beliefs' based on model 
integrable systems. Two related points emerge (i) the 
belief that the number of constants of motion is ~ L cx 
log (Af) where L is the number of sites and Af the size 
of the total Hilbert space and (ii) the so called 'rule of 
three', i.e. the belief that any many body lattice model 
in 1-dimcnsion with a fixed number of particles reveals 
its integrability only in the three particle sectoJ^, since 
the one particle sector and two particle sectors have as 
many constants of motion as the particle number (total 
energy and momentum). 

With regard to (i), we distinguish between the much 
larger Af and N of this work. By a process of block di- 



agonalizing the Hamiltonian operator into different sec- 
tors, one arrives at a direct sum representation of the 
full Hamiltonian. Each sub block is 'irreducible' in the 
sense that the space-time and internal space (parameter- 
independent) symmetries have been extracted out, and 
our considerations revolve around such sub blocks with 
a smaller and variable dimension N. Our point is that 
any such sub block must be special in the sense discussed 
here. Stitching back the irreducible blocks to reconstruct 
the full Hamiltonian matrix requires a detailed knowledge 
of the symmetries used in the first place. While possible 
in principle, we regard this process as of secondary impor- 
tance as compared to the one undertaken here, namely 
the characterization of the sub blocks themselves. 

With regard to (ii) our studies of two typical exam- 
ples give some insight into this question. Firstly, the 
Gaudin magnets show that the two and higher particle 
number sectors yield matrices that have non linear con- 
servations laws in addition to the linear ones discussed 
here. Secondly, we can study the fermionic representa- 
tion of type 1 matrices, Eq. (351. Here n = L = N, the 
sector n e = 1 is isomorphic to type 1 matrix family (JsJ) , 
while other sectors are much larger matrices of high types 
M = ( w ) — JV+lE3 i.e. they are apparently distinct and 
much more complicated integrable matrix families. With 
the hindsight of Eq. ( 35 ) it is evident that they all are 



different manifestations of the same type 1 model and 
are in some sense equivalent, but can this be formulated 
generally so that all such similar matrix integrable mod- 
els are naturally grouped together and recognized to be 
related to each other as different representations of the 
same underlying structure? Further work is needed to 
obtain clarity on these questions. 

Other open questions include: a general construction 
of type M > 3 integrable families, analytical results for 
crossings in types M > 1 and for level statistics in all 
types, the relationship between the exact solution for 
types M > 1 through a single algebraic equation and 
Bcthe's Ansatz. 
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